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A New Method Applicable to Study Solid Compounds with 


Multiple Polyhedral Structures 


Zhenlian Chen and Jun Li* 


A new direct summation method, named as polyhedron 
method, is proposed to calculate Madelung energy. This 
method calculates sums of electrostatic interactions over sets 
of neutral polyhedron unit pairs rather than conventional ion 
pairs; this gives Madelung constant in a matrix. With robustly 
rapid convergence, polyhedron method is generally applicable 
for complex compounds containing multiple polyhedral 
building-blocks and numerical polyhedral connection modes. 
The matrical analysis suggests face-sharing between octahe- 
dral pairs and edge-sharing between tetrahedral pairs can be 
electrostatically stable, against Pauling’s third rule. Further, the 


Introduction 


Coordination polyhedra are important building-blocks of com- 
plex compounds, helping understand the nature of chemical 
bonding, for example, crystal field splitting, d-electron configu- 
rations, and d-p orbital overlaps." =?! Their stacking, that is, the 
network of building-blocks, forms the framework of crystals 
and determines materials property. For a given closed packed 
anion sublattice, the cation distribution determines the con- 
nection of polyhedra. Studies of those polyhedral structures 
not only provide insights on complicated phenomena such as 
phase transition but also offer a simple yet high-throughput 
design of compounds that can be infused with different ele- 
ments to modify their physical and chemical properties.” 
For example, solid-solution designs are a major approach in 
pursuing better electrochemical performances in lithium-ion 
cathode compounds.®'® During delithiation and lithiation 
processes, the intermediate delithiated phase of cathode com- 
pounds often involves the energy competition of occupation 
and vacancy on Li sites, accompanied by the charge order of 
transition metal ions in mixing valence states. 

Although great progresses have been achieved in predicting 
crystal structures merely from knowledge of chemical compo- 
sition, there remain two challenges in materials design." 1!® 
One is the combinational complexity of searching strategy; the 
other is the energy landscape of global minimization. A closed 
packed anion sublattice with cations on interstitial sites (octa- 
hedral and/or tetrahedral) provides a starting platform to 
design new solid compounds by construction.’ Within an 
anion sublattice, only a finite number of cation distributions 
can be allowed by criteria of geometric restriction and energy 
filtering."”1 The design becomes straight: first, all possible 
polyhedral units can be generated by a combinator; then, non- 
repeating configurations of cation distributions can be filtered 


matrical calculation of Madelung energies offers a unique 
advantage to evaluate enormous configurations of cation dis- 
tributions in a given lattice in a high-throughput manner. That 
is applicable to study solid solution composites, polymor- 
phism, and defect structures, including but not limited to 
intermediate phase of delithiated cathode compounds, charge 
order or antisite defects, and extensively magnetic order. 
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and grouped by rapid energy calculations. The aim of the 
present work is to invoke Madelung energies (MEs) calcula- 
tion as an efficient filtering algorithm for global energy land- 
scape. In coupling with first-principles calculations, details of 
further structural modifications and material properties can 
be obtained on selected candidates without loss of general- 
ity. This approach has a physical root on that ME is the elec- 
trostatic energy within point-charge hypothesis, an overall 
distribution of all elements in a chemical compound. As a 
major factor to lattice cohesive energy, ME is also a basis to 
understand a variety of physical and chemical properties, 
including piezoelectricity and electrochromism, even protein 
folding." 8?" 

While supercomputing has advanced in the area of compu- 
tational chemistry, the calculation of ME is essentially the same 
as at the time Ewald Summation was proposed. It is still time 
consuming and very inefficient to calculate ionic configura- 
tions in a huge number. Madelung constant (MC), the ratio 
between the electrostatic energy due to the presence of the 
periodic images and that due to the lonely ion pair, has to be 
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Figure 1. Schematic illustration of polyhedron method. a) Crystal lattice of NaCl and neutral NaCl, NPU. Red big balls are Na ions, black small Cl. b) Illustra- 
tion of a NaClg NPU pair. Here, rij, fp, and rm are Na-Na ions vectors (from center cation of NPU-i to that of j), Na-Cl ionic bond vectors within NPU-i (from 
cation center of unit į to its anion corner n), and within NPU-j (from cation center of unit j to its anion corner m), respectively. And, o and qj are the 
charge of cations in unit i and j, respectively. NL and NM are coordination number of cations in unit i and j, respectively. c) Group Coulombic interactions 
U between NPUs, where subscript C is for cations, A for anions and the summation is over all corners of each NPU. 


calculated on each possible structure, no matter what meth- 
ods, integral transformation or direct summation methods, are 
used.24-2°! The difficulty roots in the definition of MC, which is 
a function of lattice, ion coordinates, as well as valence states 
of ions. 

Here, based on Pauling’s rules of electrostatic valence and 
taking ionic crystalline structures as a periodic piling of coor- 
dinating polyhedra," we developed a new method, named 
as polyhedron method, to calculate MC. The method counts 
Coulomb interactions between polyhedron pairs instead of 
ion pairs. That results in very rapidly robust convergence 
without inviting any mathematic trick to deal with boundary 
conditions of lattice summation and is easy to implement. 
Notably, the expression of MC adopts a matrix form directly 
correlating electrostatic potential to polyhedral structure. The 
matrix is a function of lattice and ion coordinates only, 
excluding valence states of ions. That facilitates high- 
throughput calculation of MEs for all cation distributions 
within the same anion sublattice. 


Methods 
Description of the polyhedron method 


A. Construction of Neutral Polyhedron Unit. The core imple- 
mentation of polyhedron method is the construction of neu- 
tral polyhedron unit (NPU), which is illustrated with example 
ionic crystal NaCl, shown in Figure 1a. Where, a full point 
charge (+1) is assigned to the center cation Na, while a frac- 
tion charge (—1/6) to coordinating anion Cl. The neutrality of 
the unit guarantees the network assembled by NaCl octahe- 
dral NPUs remains neutral and the construction of NPU could 
be generally and easily implemented for any complex com- 
pounds within point charge approximation. 


B. Representation of MM and ME. Figure 1c gives a set of Cou- 
lombic interactions among all ionic pairs between a NPU 


pair with parameters shown in Figure 1b, and the group 
Coulombic interaction Uu between NPU-i in center cell and 
NPU-j in L cell (L is the lattice vector given by 
L= ha+hb+hçc, a, b and c the unit cell vectors) is: 
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where ry = rj+L and . 

Because the neutrality of polyhedron unit radically removes 
net point charge oscillation in lattice summation, the summa- 
tion of Ui over lattice L is now unconditionally convergent to 
a definite physical result for any given numerical threshold, 
while conditionally convergent in original format(the specified 
condition to obtain a definite result is shown in Appendix A 
and B). A matrix form independent of the point charge can be 
obtained: 


1 
M; = UL, (2) 
f E 


No boundary condition, no special lattice cut-off, even number 
of ions arrangement along expanding edge, charge compensa- 
tion on the summation boundary or average over nanoparticle 
is employed in the polyhedron method.?*773 The conver- 
gence of each matrix element Mj is solely determined by the 
geometric combination of NPUs in the context of periodic 
arrays. Note, the matrix (Madelung Matrix, MM) is different 
from conventional definition of MC, which is referred to indi- 
vidual ion depending on the point charge assignment in com- 
pounds and is inconvenient for multicomponent complex 
crystals.” It should be noted that MM is in the unit of the 
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NPU 1 (referenced NPU) 


Figure 2. NPUs and cell expansion for NaCl calculation. a) Four NPUs in a cubic cell. b) NPUs involved with matrix elements M,, and Mu, d Three dimen- 
sional lattice expansion (nL). The OL cell is the center cell, represented by the red cube. 


reciprocal of the distance, different from unit-less MC; only cat- 
ions within unit cell enter matrix index because the calculation 
entity is the whole NPU group, not individual ions as imple- 
mented in conventional MC approach. 

For a generic compound with an array of cationic charges 
Q= (qi, dn... qn) per unit cell, ME is calculated through the 
quadratic form of MM: 


1 ' 
Eu = 5 QMQ, (3) 


where Q' is the transpose of cationic array Q and the factor 2 
removes double counting of NPU pairs. Once the matrix is cal- 
culated, ME for arbitrary cationic distribution within that struc- 
ture family can be evaluated in a high-throughput manner by 
matrix algebra of eq. (3), providing calculation efficiency for 
charge ordering determination for cations with variable 
valence states and materials of variable stoichiometry. The 
conventional MC a can be transformed from ME by «=—Em/N¢ 
XR or a=—Em/NeXa (Nç is the number of formula units in the 
unit cell), referred to bond length R or lattice constant a, 
respectively. 


C. Calculation Procedure of NaCl as Example. Here, uses con- 
ventional NaCl unit cell (with four Lu. per cubic cell) to illus- 
trate the whole procedure of Polyhedron Method. Four NPUs 
are constructed on the four Na sites, seeing in Figure 2a, giv- 
ing 4 X 4 MM. The four NPUs are equivalent, then all diagonal 
(Mi, i = 1,4) and non-diagonal (ad elements are identical, 
respectively. The two polyhedral building-blocks relating to 
distinct elements Mı and Miz are shown in Figure 2b. The 
interactions contributing from the center cube to Mi, Mo, 
are 


Ucc=0, 


GEET D RE 
Uca Uca 


1 1, 4 
Uu — X4 — +—— $X6 
M 36 in =} 


where, d (2.8201 A) is the Na-Cl ionic bond length. Then, 
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Three dimensional cell expanding along unit cell vectors 
from the center cell is shown in Figure 2c and the summation 
is over nL cells, (n + 1)X(n + 1)X(n + 1) of cubes at summation 
length nL (n=0,1,...; L is a symbol representation of the 
length of unit cell vectors). The evolutions of calculated values 


Table 1. MM elements (unit: A) as functions of summation length nL. 


Shell Nusty mnt mnt a 
OL 1 —0.5125 —0.0433 0.9059 
1L 3°-1=26 SC —0.1932 —0.1343 

TC —0.7057 —0.1777 1.7466 
2L 53°—3°=98 SC 0.0001 —0.0002 

TC —0.7056 —0.1779 1.7476 


SC and TC represent contributions from cubes in the nL shell and all 
cubes up to nL cell, respectively. The « is the converted MC at the sum- 
mation length nL. Nypus is the number of NPUs in the nL shell. 
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Table 2. MEs Eu (unit: eV/f.u) as a function of summation length nL for cubic NaCl and CsCl, Wurzite ZnO, undistorted and distorted perovskites SrTiO3 


and LaMnO3 and rhombohedral cathode LiCoO3. 


oL 1L 2L 3L CL E am al Rora 
NaCl —4.6256 —8.9183 —8.9234 —8.9234 2L 1.7476 1.747568"! 2.8201 
CsCl —2.789 -7.1273 —7.1085 —7.1087 EI 1.7627 1.76268?" 3.5706 
ZnO —24.7536 —47.6913 —47.9881 —47.9953 4L —47.9940 5.9872 5.994138) 1.7964 
StTiO3 — 104.3587 — 182.5812 — 182.5718 —182.5716 3L 49.5099 49.5122?! 3.905 
LaMnO; —107.757 — 166.6544 — 166.6967 —166.6996 4L —166.7007 22.0195 1.9021 
LiCoO, —57.9281 —108.1180 —111.0780 —111.0850 3L 38.1554 4.9461 


CL is the convergence length of all matrix elements reaching the converging threshold (TE — E; 
of ME and MC calculated by the polyhedron method, respectively. « 


ref 


EI TT Lu) < 1 fe 7. E and gP™ is the converged values 


is the MC value reported in literature. R or a are values used in the reference MC. 


of matrix elements and conventional MC at summation length 
nL are listed in Table 1. At summation length 2L, the values of 
MM reach the preset converging threshold (cf. Table 2) and 
give the conventional MC 1.7476 by the conversion 


1 1 
o qfux2.8201 =- g QX MxQ'X2.8201 


=0.7056. —0.1779 —0.1779 =0.1779 
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—0.1779 —0.7056 —0.1779 —0.1779 1 
=-—0.3525X(1,1,1,1)X x 
0.1779 —0.1779 —0.7056 —0.1779 1 


=0:1779. =0.1779 
=1.7476. 


—0.1779 —0.7056 


Notably, Table 1 shows the contribution from the 2L shell 
has become small in order of 1074 A"! for both diagonal ele- 
ment and non-diagonal element. Here, nL shell means all 
cubes locate between nL and (n—1)L length. At summation 
length 2L, that is, 4d, the 2L cell consists of 5 X 5 X 5 X 
8 = 1000 ions, in contrast to 8L, that is, 16d or 27000 ions 
required in Bakers’ scheme (cf. Table 3 in Ref. 29), which was 
the shortest length to give the value 1.7476 by direct MC cal- 
culation in literature to our best knowledge. 


Solid-solution configurations construction 


We take spinel LiMnz0, as example to illustrate the determination 
of the cation distributions using ME as a high-throughput filter. 
The flow chart is shown in Figure 3. In spinel LiMnzO,, Mn exhibits 
variable valence states, Mn?* and Mnî*. In Fd3m unit cell, there 
are 16 Mn ions on octahedral sites and 8 Li ions on tetrahedral 
sites, generating a 24 X 24 MM. When one picks up 8 octahedral 
sites for Mn? ions among the 16 octahedral sites, the combinator 
gives C$, 12870 configurations. Each can be represented by a cati- 
onic array Q (e.g. [1,1,1,1,1,1,1,1,3,3,3,3,3,3,3,3,4,4,4,4.4.4.4,4], where 
the first eight for Li ions and the other 16 for Mn ions). Then by 
comparing MEs calculated for the 12,870 set by eq. (3), 25 distinct 
configurations are obtained in a high throughput manner. Thus, 
this method can be generally applied to other solid solutions of 
elements even in identical nominal valence state or antiferroma- 
gentic order, only if elements or spin orientations are assigned to 
different pseudo valence states, which are not necessary to have a 
chemical meaning but will help identify or classify geometric pat- 
terns of cation distributions. Detailed chemical bonding effects can 
be captured by following-up first-principles calculations. 


Applications 
General applicability of the polyhedron method 


The neutrality of polyhedron unit ensures the polyhedron method 
a robust approach to MC calculations, regardless of Bravais lattice 
type, different from other reported direct summation methods, 
which are most applicable to cubic compounds.””*"! A systematic 
calculations of ME for various compounds, such as often-studied 
cubic NaCl and CsCl, Wurzite ZnO, perovskites SrTiO3 (undistorted, 
cubic), and LaMnO3 (distorted, orthogonal) and rhombohedral 
cathode LiCoO;, all reach flatting values within numeric precision 
0.1% at a summation length not longer than 2L. Digital precision is 
slightly structure dependent but ME always converges monotoni- 
cally, seen in Table 2. The accuracy of all converged MEs at the con- 
verging threshold is better than chemical accuracy. From the point 
of view of a high efficient and high throughput energy filtering, 
we do not pursue further accuracy and precision in the calculation 
of MC for structure prediction beyond point-charge hypothesis. 


Physical insight of polyhedral connection 


What MM describes is the Coulomb interaction between polyhe- 
dral building-blocks, analog to Pauling's structural characteristics of 
coordination compounds, especially relating to his third rule about 
polyhedral combination. Figure 4 shows the stacking of all octahe- 
dral and tetrahedral interstitial sites within cubic close-packing (AB 
stacking) and hexagonal close-packing (ABC stacking) oxygen sub- 
lattices. One cation-oxygen octahedra layer sandwiched between 
two cation-oxygen tetrahedra layers. There could be nine connec- 
tion modes: corner-, edge-, and face-sharing between tetrahedron 
and/or octahedron in AB close-packing oxides; whereas six modes: 
corner- and edge-sharing between octahedron-octahedron pair 
and tetrahedron-tetrahedron pair, and corner-, face-sharing 
between octahedron-tetrahedron pair in ABC stacking. The MM 
elements for the modes (values are given in Table 3) can be 


Table 3. MM elements M; for NPU pair (i, j). 


Tetrahedron- Octahedron- 
Octahedron-octahedron tetrahedron tetrahedron 
Pairs 
sharing corner edge face corner edge face corner edge face 
WON 0.91 1.00 0.71 058 0.71 0.5 0.29 083 06 0.43 


MÂ’) 0.17 —0.19 —0.24 —0.24 —0.32 —0.41 0.42 —0.24 —0.31 —0.13 


Here, rij is the distance between the centers of the pair. 
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Figure 3. The flow chart illustrating construction of solid solution configu- 
rations using ME as a high-throughput filter algorithm. [Color figure can be 
viewed in the online issue, which is available at wileyonlinelibrary.com.] 


calculated from a hexagonal supercell of six layers with lattice 
parameters a = 3d9 and c = 4vV3ao, where do = 0.707 Å. When 
the six stacking of oxygen layers is ABABAB, the cation-oxygen 
polyhedral sequences are AAAAAA and ABBAABBAABBA for octa- 
hedra and tetrahedra, respectively; for the ABCABC stacking, the 
sites are ABCABC and ABCABCABCABC for octahedra and tetrahe- 
dra respectively. For the six polyhedral connection modes in ABC 
stacking, the distances between cation centers of polyhedra are 
the same as those in AB stacking, except for octahedron- 
octahedron corner-sharing, of which the distance is 1.00 Å and 
0.91 Å in ABC and AB stacking, respectively. 

One finds that with tetrahedron, face-sharing is often more 
electrostatically instable. Especially, face-sharing between 
tetrahedron-tetrahedron pair even gives a positive value, confirm- 
ing Pauling’s speculation that face-sharing intend to decrease the 
stability of crystals with low coordination number However, 
between octahedron-octahedron pair, face-sharing gives value 
extremely close to edge-sharing and negatively larger than 
corner-sharing; among all polyhedron pairs, edge-sharing gives 
values negatively larger than corner-sharing. That indicates Paul- 
ing's third rule partially incorrect and supports that edge-sharing 
between tetrahedra can stabilize structures of /-BeO, dumortier- 
ite Al,BO3(SiIO,)303, LisFeO, and LigCoO,, and face-sharing 
between octahedra can stabilize 2H-BaMnO3.2? 3" The physical 
root is that the collective Coulomb attraction between anions of 
NPU-j and cation in NPU-i surpasses cation-cation Coulomb repul- 
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Figure 4. Connection modes between polyhedral building-blocks in AB and 
ABC stacking of oxygen sublattices. Black small, Green large and purple 
small balls represent oxygen, octahedral (octa) sites, and tetrahedral (tetra) 
sites, respectively. Labels C, E, and F represent corner-, edge-, and face- 
sharing, respectively. [Color figure can be viewed in the online issue, which 
is available at wileyonlinelibrary.com.] 


Table 4. Energy terms and magnetic moments for the four high symme- 
try charge-order configurations. 


Energy Magnetic moments 
Eo E 5 Mn2+ Moi" oää oa: 
—247.72 —47.321 1/2 3.63~3.66 —0.09~-0.11 
Imma —246.68 —47.515 1/6 3.89 3.15 0 —0.08 
P4322 —246.72 -—47.452 3.9 3.17 —0.05 004 
Cmm2 —246.96 —47.507 3.9 3.15, 3.18 ~-0.01 ~-0.08 


Em and E are MEs and our first-principles GGA + U total energies by 
VASP structural relaxations, respectively.2***! For the reference Fd3m 
configuration with all Mall", ME is —244.58 eV/f.u. Magnetic 
moments are calculated in ferromagnetic order, representing charge 
disproportion on oxygen sites. 


sion for most connection modes except for face-sharing between 
tetrahedron pair, but it was ignored in Pauling’s speculation. 


Charge order determination in coupling with first-principles 
calculations 


Now, we use spinel LiMn 20, again to show that in coupling with 
first-principles calculations, polyhedron method can provide 
more deep insights on the physical and chemical properties for 
complex compounds. There are four configurations with high 
symmetry (by ISOTROPY analysis?!) of P4m2, P4322, Imma, and 
Cmm2, among the twenty-five distinct configurations according 
to MEs. In Fd3m lattice, all MnO, octahedra are undistinguishable 
then the major contributions, intra-NPU interactions of all MnO. 
units, are the same. The ME for Mn**-Mn** charge-order is 
always lower than the referenced Mn®!** configuration (cf. 
Table 4). This is due to the quadratic ruling in Coulomb energy 
[referring to eq. (3)], cf., sum of square (37 + 47) vs. square of aver- 
age (2 X 3.5’). This simple rule indicates charge orders always 
electrostatically stabilize a compound with mixing valences if dif- 
ferences in NPUs surpass other electronic repulsion effects. 

In Mn?*-Mn** charge-order configurations, charge dispropor- 
tion forms on oxygen ions because each oxygen ion is shared by 
one UO, tetrahedron with three MnOg octahedra. In P4m2 config- 
uration, three Mn?* or three Mn** ions share one oxygen ion, 
then the valence disproportion is 1/2 and the nominal valence of 
oxygen is -(2-1/2) or -(2 + 1/2); whereas, two Mn?*, one Mn** or 
two Mn**, one Mn?* share one oxygen ion in the other three 
configurations, then the valence disproportion of oxygen 
becomes smaller to 1/6, seen in Figure 5. In P4m2 configuration, 
the valence disproportion is so strong that it is leaning to an aver- 
age valence state after first-principles structural relaxation; 
whereas the other three configurations all favor clear Mn? *-Mn** 
charge-orders. These features are revealed well in magnetic 
moments patterns, seeing Table 4. Total energies of the three 
charge-order configurations are again lower than that of P4m2 
configuration. The atomistic difference among the three Mn?™- 
Mr?" charge-order configurations lays in the pattern of Mn’! " 
O,/Mn"! + Os combinations, that is, the nearest Mn-Mn connec- 
tions, cf. Figure 5, and gives rise to small difference in MEs as well 
as in total energies. The proximity of energy suggests various 
charge- order configurations may coexist in LiMn 0,; the charge- 
order may be local rather than global as demonstrated in the long 
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Figure 5. Landscape of MEs and MnOg octahedral combination sketched 
for charge-order configurations, where green and black small balls repre- 
sent lithium and oxygen ions, respectively. [Color figure can be viewed in 
the online issue, which is available at wileyonlinelibrary.com.] 


range partially charge-order in Fddd superlattice with 504 ions 
each cell.4! That is rendering to further experimental 
verification. 


Conclusions 


Based on the construction of NPU, a new direct summation 
method is proposed to calculate MEs for complex crystal com- 
pounds. Due to screening effects of NPUs, the convergence of 
lattice summation becomes rapidly robust and generally appli- 
cable. The MC is expressed in a matrix form and each matrix 
element accounts for the group interaction between NPU pair. 
The matrix quantitatively characterizes connection modes 
among polyhedral building-blocks. That justifies Pauling’s third 
rule and validates face-sharing between octahedra and edge- 
sharing between tetrahedra could be electrostatically stable. 
The matrix form also benefits comprehensive and systematic 
ME calculations for solid materials with variable stoichiometry 
and mixing valence states. It offers a new strategy to study 
atomistic details of solid solution and polymorphisms (such as 
A2BSiO, silicates with cations A, B, and Si on tetrahedral sites 
in AB stacking oxygen lattice), charge order, or even spin- 
polarized systems, for example, anti-ferromagnetic order. This 
may guide the searching of charge-pairing mechanism in very 
challenging solid compounds and provide a new strategy to 
study complicated processes for compounds with variable 
valence metal ions including intercalation cathode of lithium- 
ion batteries. 


APPENDIX A 


Mathematic Proof of Unconditional 
Convergence of Lattice Summation of Uj 
Due to NPU construction, net point-charge interactions are 


screened off between NPUs, and eq. (1) can be rewritten as 
sum of interactions between dipolar pairs: 
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The first three terms in eq. (A2) are dipolar contributions and 
the last term includes all other higher multipolar terms with 
smaller contributions. The dipole terms vanish after merging: 


C m=1,N- SE i 
1 (A3) 
gr D Val m—¥n)=0 
Then, 
L eri 


1 1 
3 ro( 2 )=0( 3): (A4) 
Io éi Ti 


That means Ui attenuates to the third order or higher order of 
reciprocal of distance between the centers. Therefore, the lat- 
tice summation of LE should converge monotonically, scaling 
as 1/rj because the number of polyhedra increases with cell 
expansion scaling as ry In other words, for any given con- 
verging threshold A, one can find a suitable summation length 
nL, satisfying WEI EU Tip < A. That confirms the con- 
vergence of the polyhedron method is mathematically robust, 
in contrast to the conditional convergence of classical MC 
schemes. If the symmetry of the polyhedra is high, higher mul- 
tipolar terms such as quadrapolar term may also be counterac- 
tive so as to achieve much faster convergence than dipole 
convergence. 


APPENDIX B 


Independence of Convergences on Cell 
Expansion and Lattice Shape 


The neutrality of the polyhedron unit is a strong physical 
imposition on the group summation that drives the conver- 
gence of the polyhedron method, insensitive to the boundary 
condition, that is, the scheme of cell expansion in sphere or 
cube and lattice shape in cubic or triclinic. 

Figure B1 shows four cases using NaCl as example for sphere 
and cube cell expansion with primitive cell (PC) and crystal 
cell (CC), all cases reach convergence within 5L with error 
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Figure B1. The calculation of ME for NaCl. a) Unit cell of PC and CC. b) The cations and polyhedra on boundary in 1L summation using sphere expansion 
with PC. d The convergence using sphere and cube expansion with PC and CC. [Color figure can be viewed in the online issue, which is available at 


wileyonlinelibrary.com.] 


smaller than 10 +. This is in contrast to the conventional direct 
summation method, which is divergent for sphere lattice 
expansion in the calculation of MC for NaCl without correction 
term.2°! In conventional methods, whether or not the ions, 
cations or anions, are on the boundary is determined by its 
pair distance to the refer Jon Hl On the contrast, in polyhe- 
dron method, the determination condition is the pair distance 
between the center cations of the polyhedron unit and the 
refer unit. In fact, coordinated anions are on the boundary. For 
example, shown in Figure B1 b, the Cl ions coordinated to the 
12 cations (forming icosahedron) seen in the upper atomic 
model in Figure B1 b are on boundary of the 1L summation 
with PC using sphere expansion. 

While the value of ME using sphere expansion with rhombo- 
hedral PC deviates far away from the convergent value in sum- 
mation length of 1L, it converges rapidly to the same definite 
physical result after 5L. That deviation roots in the group of 
cations, which is determined by the lattice constant in our cal- 
culations. In the sphere expansion with PC, the radius of 1L 
length is of 3.9882 A and only nearest cation neighbors are 
included; while at least second nearest cation neighbors are 
counted in for the other three cases. 

Appendix A has shown that polyhedron method will obtain a 
definite result,. which is independent on the boundary scheme 
of cell expansion and lattice shape. However, this does not 
exclude the numeric deviation before convergence reaches. 
The numeric deviation of the polyhedron method relies on the 
symmetry of the polyhedron units and the lattice. The latter 
one determines the networking of the polyhedra on boundary, 
kind of similar to other direct summation methods. This devia- 
tion is different in nature from the convergent difficulty in 
other direct summation methods, in which the convergent dif- 
ficulty or errors mainly come from the correction terms on 


boundary. Recent works showed such difficulties also exist in 
the well-developed Ewald summation.4* 
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A new method is proposed to study electrostatic interactions among polyhedral building-blocks in 
complex crystals within point charge approximation. By counting group interaction between neutral 
polyhedron unit pairs, the Madelung constant is obtained in matrix form, revealing the geometric cor- 
relations among interstitial sites in lattice. Distinct cationic distributions can be classified according to 
Madelung energies by high-throughput matrix algebra. 


